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UNSTEADY  TRANSONIC  FLOW  COMPUTATIONS 


A.  R.  Seebass,  N.  J.  Yu,  and  K-Y.  Fung 

Department  of  Aerospace  and  Mechanical  Engineering 
University  of  Arizona 
Tucson,  Arizona  85721 


SUMMARY 


We  Investigate  the  effects  of  unsteady  modes  of  motion  on 


Invlscld  small  perturbation  approximation.  The  study  Is  a numerl 


cal  one  and  draws  upon  the  alternating-direction  Implicit  procedure 


developed  for  such  calculations  by  Ballhaus  and  his  coworkers  at 


the  NASA  Ames  Research  Center.  Our  numerical  algorithm  treats 


shock  waves  as  moving  discontinuities.  Results  of  nonlinear  and 


time- linearized  calculations  of  the  transonic  flow  past  an  NACA 


64A006  airfoil  experiencing  harmonic  motions  In  several  of  Its 


1.  INTRODUCTION 


In  unsteady  transonic  flows,  relatively  small  periodic  changes  In  the 
boundary  conditions  can  lead  to  substantial  changes  In  the  loads  and  moments 
with  marked  phase  lags.  These  are  of  major  concern  In  the  aerod3mamlc  de- 
sign of  aircraft  that  operate  In  the  transonic  regime.  A short,  but 
timely,  review  of  various  aspects  of  unsteady  transonic  flow  may  be  found 
In  Reference  1.  Of  particular  concern  are  aero-elastlc  behavior,  and  flutter 
and  buffet  boundaries.  Here  the  unsteady  perturbations  may  sometimes  be 
small  enough  that  linearization  about  a nonlinear  steady  flow,  as  suggested 
by  Landahl  (2)  long  ago.  Is  possible. 

In  such  flows  the  behavior  of  the  boundary  layer,  especially  as  It  Is 
affected  by  the  pressure  rise  caused  by  any  shock  waves  In  the  flow.  Is 
clearly  of  major  Importance.  Additionally,  In  the  neighborhood  of  the 
leading  edge  the  flow  perturbations  are  large;  consequently,  highly  accurate 
Invlscld  results  require  the  use  of  the  full  potential  equation.  It  seems 
likely  that  eventually  the  computational  algorithms  used  In  routine  studies 
of  unsteady  transonic  flows  will  use  the  Reynolds  averaged  Navler-Stokes 
equations  now  used  In  research  studies.  However,  such  algorithms  (3)  cur- 
rently require  substantial  computer  time  and  are  too  Inefficient  for  explora- 
tory studies  such  as  this  one.  The  ability  of  these  algorithms  to  model 
complex  unsteady  transonic  flow  phenomena,  such  as  buffet,  has  recently  been 
demonstrated  (4) . 

An  important  consideration  In  constructing  an  algorithm  for  unsteady 
transonic  flows  Is  the  treatment  of  moving  shock  waves.  The  experimental 
observations  of  Tljdeman  (5-7)  Indicate  that  even  for  simple  airfoil  motions 
shock  wave  motions  can  be  complicated,  and  that  they  can  strongly  affect 
aerod3mamlc  force  and  moment  variations.  Time- linearized  methods,  l.e.. 


methods  that  assume  the  unsteady  perturbations  are  small  compared  to  the 
basic  steady  disturbance,  have  not  usually  considered  shock  motions  (8,  9), 
although  they  can  be  modified  to  do  so  for  small  shock  excursions  (10) . 

Time- Integration  methods  (11-18)  treat  shock  waves  by  "capturing"  them,  a 
procedure  that  can  present  a number  of  difficulties. 

Unsteady  experiments  (5-7),  analysis  (10)  and  numerical  studies  (10)  all 
Indicate  that  the  amplitude  of  the  shock  wave  motions  Increases  Inversely 
with  reduced  frequency.  Thus  some  of  the  most  Important  effects  occur  with 
low-frequency  motions.  This  Is  not  surprising;  nonlinear  behavior  Is  sup- 
pressed at  higher  frequencies,  with  the  small  perturbation  equation  becoming 
linear  for  frequencies  higher  than  the  two-thirds  power  of  the  airfoil's 
thlckness-to-chord  ratio.  Explicit  finite-difference  schemes  are  not  effi- 
cient when  applied  to  lo^^frequency  cases  because  the  stability  restriction 
on  the  time  step  Is  substantially  more  severe  than  that  required  for  accuracy. 
As  a result,  efficient  seml-impllclt  methods  (13)  and  even  more  efficient 
fully  implicit  methods  (11,  12,  17,  18)  have  been  developed.  Caradonna  and 
Isom  (17)  use  an  Iterative  Imp^^lt  procedure,  l.e.,  the  nonlinear  Implicit 
finite-difference  equations  must  be  solved  iteratively  at  a given  time  level. 
In  an  earlier,  unpublished,  study  we  also  used  such  a procedure.  Ballhaus 
and  Steger  (11)  and  Beam  and  Warming  (18)  constructed  more  efficient  algor- 
ithms that  solve  the  nonlinear  equations  directly  by  the  solution  of  simple 
matrix  equations  generated  by  an  alternating-direction  implicit  (ADI)  proce- 
dure. This  method  has  proven  to  be  so  efficient  that  It  Is  now  used  as  an 
alternative  to  successive  line  over-relaxation  (SLOR)  for  steady  flow  calcu- 
lations (Reference  19  and  Yu  and  Seebass,  unpublished). 

As  mentioned  above,  these  Implicit  schemes  "capture"  shock  waves,  l.e., 
shock  waves  evolve  automatically  as  part  of  the  numerical  solution.  Shock 
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capturlng  produces  shock  profiles  chat  are  distorted  In  a manner  than  de- 
pends on  Che  truncation  errors  In  the  finite-difference  scheme.  The  use  of 
mixed-difference  schemes  (11,  18)  can  Improve  the  situation  for  cases  In 
which  the  flow  changes  from  supersonic  to  subsonic  across  the  shock.  How- 
ever, when  this  condition  Is  not  satisfied  the  differencing  cannot  be 
switched  across  the  shock  and  shock  resolution  Is  poor.  In  any  case,  shock 
capturing  requires  spatial  grid  spaclngs.  In  regions  where  shock  waves  are 
anticipated,  that  are  sufficiently  small  to  resolve  the  shock  waves.  The 
grid  spacing  required  to  do  this  Is  usually  much  smaller  than  that  required 
to  resolve  flow  variable  gradients  In  most  of  the  rest  of  the  flow  field. 
Shock  fitting  removes  the  large  gradients  from  the  finite  difference  solution 
and  permits  equivalent  flow  field  resolution  with  fewer  grid  points,  both  In 
space  and  time  (20,  21).  If  shock  waves  are  not  treated  as  discontinuities, 
but  are  to  be  captured  correctly,  the  difference  equations  must  be  solved  In 
conservation  form.  This  Imposes  an  additional  constraint  on  the  construction 
of  finite-difference  schemes  that  can  be  difficult  to  satisfy. 

A need  for  shock  fitting  also  arises  In  computing  tlme-llnearlzed  solu- 
tions for  very  small  unsteady  perturbations.  Tlme-llnearlzed  solutions  for 
Indlclal  motions  can  be  used  to  determine  force  and  moment  coefficient  vari- 
ations at  various  reduced  frequencies,  obviating  the  need  for  a numerical 
solution  at  each  reduced  frequency  (see,  e.g..  Reference  22). 

Traci,  et  al.,  (8,  23,  24)  have  developed  relaxation  methods  for  solving 
the  resulting  tlme-llnearlzed  equations  of  motion  for  harmonic  disturbances. 
Less  complete,  but  comparable,  studies  have  been  made  by  Weatherlll  et  al. 
(9);  these  derive  from  an  earlier  study  by  Ehlers  (25).  In  both  of  these 
studies  shock  motions,  which  contribute  substantially  to  the  time-varying 
loads  and  moments,  are  neglected.  Difficulties  also  arise  In  the  convergence 
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o£  the  iterative  numerical  scheme.  Unsteady  small  amplitude  motions  have 
shock  wave  excursions  that  are  the  order  of  the  amplitude  of  the  motion 
divided  by  the  reduced  frequency  of  the  motion.  Consequently,  these  shock 
motions  dominate  other  low  frequency  contributions  to  Che  lift  and  moment 
coefficients.  Such  tlme-llnearlzed  shock  motions  can  be  computed  In  a 
rational  way,  but  the  accuracy  of  the  results  depends  critically  on  an 
accurate  resolution  of  the  steady  flow  field  In  the  vicinity  of  the  shock 
wave  (10);  this  Is  best  accomplished  by  shock  fitting. 

This  paper  briefly  reviews  the  numerical  procedures  we  have  developed 
for  computing  nonlinear  and  tlme-llnearlzed  small  perturbation  unsteady 
transonic  flows.  We  use  an  ADI  scheme  and  treat  shock  waves  as  discontinu- 
ities In  Che  flow.  Calculations  of  Che  transonic  flow  past  an  NACA  64A006 
airfoil  experiencing  harmonic  or  Indlclal  pitching  and  flap  oscillations  are 
discussed. 

2.  FORMULATION 

We  write  the  unsteady  small  disturbance  equation  for  low  frequency 
transonic  flows  in  the  commonly  used  form 

+ (1  - - (Y  + V 

The  spatial  coordinates,  the  time,  and  the  velocity  potential  In  (1)  have 
been  non-dlmenslonallzed  by  the  chord,  the  reciprocal  of  the  angular  fre- 
quency, and  the  free  stream  velocity  times  the  chord,  respectively.  Ocher, 
perhaps  more  useful  and  suitable,  forms  are  given  In  References  21  and  26. 

This  equation  results  from  a systematic  expansion  of  the  velocity  potential 

2/3 

In  Che  thickness  ratio  T and  applies  for  reduced  frequencies  R > 0(t  ) 
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where  K * uic/U,  l.e. , Che  angular  frequency  nniltlplled  by  the  time  it 
takes  Che  flow  to  traverse  Che  airfoil  chord.  Lin,  Relsner  and  Tslen  (27) 
showed  chat,  with  restriction  to  small  perturbations  throughout  the  flow, 
this  is  Che  only  nonlinear  equation  that  arises.  For  moderate  frequencies 
the  equation 

^ - 2K(j»  ^ + {1  - - (y  + 1)M^[((>  + ^ ^ K(j)  j}d)  +6  - 0 

tt  ^xt  00  o»'-'*'x  Y + 1 t ^xx  ^yy 

Is  frequently  used,  with  or  without  the  term,  and  may  provide  results 

that  apply  at  higher  frequencies  than  those  obtained  from  (1) . 

The  boundary  condition  on  the  body  Cakes  the  simple  form 


4 (x.o,t)  - t[y«  +7  + ky;|)], 

y X X X u 


1 

- Jix 


(2) 


where  Y(x,c),  the  Instantaneous  body  shape,  has  been  decomposed  Into  a 

steady  part,  Y*,  and  an  unsteady  part,  y“.  Here  6 is  the  amplitude  of  ! 

2/3  ' 

the  unsteady  oscillation.  Because  K - 0(t  ),  the  last  term  in  (2)  is 

dropped  unless  Y^*  = 0 or  is  small.  For  this  reason,  the  Clme-llnearlzed 

X I 

perturbation  velocity  potential  for  plunging  motions  (Y^  - 0)  Is  Just  K 

times  that  for  Che  analogous  pitching  motion,  where  y'*(x,c)  ■ (x  - x )sln  t. 

o 

Numerical  studies  conducted  by  Magnus  (IS)  show  chat  erroneous  boundary 
data  on  a finite  domain  can  lead  to  significant  errors.  The  low  frequency 
approximation  Implies  that  any  changes  In  Che  circulation  are  communicated 
Instantly  downstream  to  Infinity.  Consequently.  *'he  simplest  boundary  con- 
ditions are  ^ downstream  boundary  and  ^ ■ 0 on  the  other 

boundaries.  Ballhaus  and  Goorjlan  (12)  used  these  boundaxy  conditions  in 
their  study  and  obtained  satisfactory  results.  The  validity  of  such  far- 

1 


field  boundary  conditions  can  only  be  justified  by  numerical  experiments; 

i.e.,  near  the  boundary  the  disturbance  quantities  it>  > and  <|)  , must  be 

X y 

much  smaller  than  the  values  at  the  airfoil  surfaces.  For  the  lifting  case, 

<p  depends  on  the  Instantaneous  circulation,  F.  This  dependence  can  be 
derived  theoretically  by  assuming  that  in  the  far  field  all  the  perturbations 
are  small  compared  to  the  basic  steady  state  (see,  e.g..  Reference  26).  Here 
we  use  a stretched  coordinate  system  that  maps  the  doubly  infinite  domain 
into  U1  i 1.  |n|  ^ 1,  and  set  ♦jj  " ® downstream  boundary  C ■ 1 

and  0 elsewhere  on  the  boundary  of  this  domain.  As  a numerical  test 

for  this  procedure  we  have  computed  the  steady  state  circulation  about  an 
NACA  6AA006  airfoil  for  various  flap  deflection  angles,  using  the  ADI  method 
with  appropriate  far'-fleld  values  of  i^,  corrected  for  the  usual  steady  state 
circulation  contribution.  These  results  have  been  compared  with  the  results 
obtained  by  the  ADI  calculations  with  the  boundary  conditions  employed  here 
for  an  unsteady  flap  deflection  to  the  correct  angle.  These  results  are 
identical  within  the  accuracy  with  which  we  have  computed  the  solutions. 

Any  shock  wave  that  exists  in  the  flow  field  must  satisfy  the  jump  rela- 
tion derived  from  the  conservative  form  of  the  governing  equation  (1),  namely 


-2KM^^  I^(dx/dt)  - (1  - - (Y  + 1)M^*  + (♦  1^  - 0 (3) 

X 8 X X Y 


together  with  the  condition  derived  from  the  assumption  of  irrotatlonality. 


(dy/dx)  - -I A 1/1 A I. 

s X y 


Here  A refers  to  the  mean  value  of  A evaluated  on  each  side  of  the 
X X 

discontinuity,  and  Indicates  the  jump  in  (p^  across  the  discontinuity; 

the  subscript  "s"  denotes  the  quantity  evaluated  at  the  shock  surface. 


■F=^ 


’1 
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The  pressure  coefficient,  defined  so  chat  It  vanishes  at  sonic  condi- 
tions, takes  the  form 


- 1 


C - -2 


(Y  + Dm; 


(5) 


In  Che  small  disturbance  approximation,  Che  Kutta  condition  Is  Imposed  by 
requiring  that  be  continuous  at  y ■ 0 for  x > 1/2. 


2.1  Tlme-Llnearlzed  Equations 

We  now  assume  that  the  unsteady  disturbances,  characterized  by  6,  are 
small  enough  chat  we  may  write 


iKx.y.t)  - (^•(x,y)  + 6i|;(x,y,t)  + o(6) 


and  neglect  higher-order  terms  In  5.  The  restriction  Imposed  on  6 for 

this  to  be  true  will  depend  on  Che  other  parameters  of  Che  problem,  viz., 

2 2 2/'? 

K 5 (1  - M^)/[(Y  + and  K.  This  gives 


(6) 


{1  - - (Y  + - 0.  ♦•(x.O)  - TY»' (x), 

A XX  yy  y 


(7) 


-2KMii|»  + {[1  - - (y  + l)Mf<>'*]4»  } + <l>  - 0, 

• xt  ' • X X X yy 


(8) 


K»  (x,o,t)  - y"  + ky;, 

y Xu 


The  solution  Co  (7)  must  satisfy  Che  steady  version  of  the  shock  relations 
(4)  and  (5).  The  shock  relations  for  (8)  are  discussed  In  Section  2.2. 

We  avoid  writing 


L 


i^(x,y,t)  - Re{^ (x,y)e^“^}  (9) 


as  this  restricts  the  study  to  harmonic  motions.  Because  indlclal  motions 
can  be  superimposed  to  obtain  the  results  for  any  frequency,  they  seem  more 
Important.  Equation  (9)  results  in  an  equation  for  a complex-valued  ij) 
which  may  be  solved  by  line  relaxation.  Our  experience  with  unsteady  ADI 
techniques  has  been  that  they  are  at  least  as  effective  as  line  relaxation 
for  problems  of  this  type,  and  hence  there  Is  no  advantage  to  the  decomposi- 
tion (9). 

The  numerical  algorithm  developed  In  Reference  21  and  described  briefly 
In  Section  3 can  be  used  to  solve  the  basic  equation  (1)  subject  to  the 
boundary  conditions  (2),  the  shock  conditions,  (3)  and  (4),  the  far-fleld 
boundary  conditions  and  the  Kutta  condition.  Steady  state  solutions, 
<P°(x,y),  may  be  obtained  rapidly  by  subjecting  a basic  steady  state,  such 
as  undisturbed  flow,  to  rapidly  changing  boundary  data  until  a new  steady 
configuration  Is  prescribed.  This,  then,  determines  the  steady  state  result 
for  (7)  needed  to  solve  (8). 


: 2.2  Shock  Fitting 

The  basic  algorithm  for  shock  fitting  In  mixed  flows  was  developed  In  a 
previous  study  of  steady  transonic  flows  (20).  A different  approach  to  shock 
fitting  has  also  been  used  by  Hafez  and  Cheng  (28)  In  their  study  of  steady 
transonic  flow  problems.  Their  procedure  essentially  replaced  the  shock- 

^ point  operator  of  Murman  (29)  by  an  analogous  difference  statement  derived 

i 

1 

I from  the  shock  jump  conditions.  Subsequently,  the  velocity  potential  on  each 

I side  of  the  shock  wave  Is  extrapolated  to  locate  the  shock  wave. 
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ITo  understand  Che  shock-fitting  procedure  for  unsteady  transonic  flow 

calculations  It  Is  necessary  to  recall  how  shock  waves  form  In  an  unsteady 
field.  Shock  waves  are  generated  when  the  local  flow  becomes  supersonic 
and  compressive.  While  the  Initial  shock  formation  may  not  be  predicted 
exactly  by  the  numerical  solution  when  shock  fitting  Is  used  In  Che  early 
stages  of  shock  wave  formation.  It  eliminates  spurious  oscillations  In  the 
numerical  solution  and  does  provide  Che  correct  development  of  Che  shock 
wave  In  later  stages  of  the  calculations  (30).  The  criteria  that  we  set 
for  the  Initial  shock  formation  Is  that  the  local  flow  become  sonic  (relative 
Co  the  airfoil)  and  compressive.  In  Che  body-fixed  coordinate  system,  a 
shock  wave  can  exist  both  In  the  usual  supersonic-supersonic  and  supersonic- 
subsonic  transitions,  but  also  In  a purely  subsonic  flow  field,  sometimes 
referred  Co  as  a "subsonic-subsonic"  shock.  In  any  case,  the  flow  ahead  of 
the  shock  relative  to  a coordinate  system  fixed  on  the  shock  must  be  always 
supersonic.  Consequently,  the  correct  judgment  for  the  existence  of  a shock 
wave  In  the  unsteady  field  Is  to  evaluate  Che  local  flow  velocity  ahead  of  a 
prospective  shock  with  respect  to  the  coordinate  system  fixed  on  It;  l.e.. 

If  Che  local  flow  Is  supersonic  a shock  may  exist.  If  the  local  flow  becomes 
sonic  Che  shock  strength  diminishes,  and  If  It  Is  subsonic  a shock  cannot 
exist. 

Any  shock  wave  that  exists  In  the  flow  field  must  satisfy  the  jump  rela- 
tions (3)  and  (4).  In  two-dimensional  small  perturbation  transonic  flows 
Che  shock  waves  that  usually  occur  are  nearly  normal  to  Che  flow  direction. 
While  It  Is  not  necessary  to  do  so.  In  the  numerical  calculations  reported 
here  we  have  assumed  Chat  If  the  basic  steady  flow  has  a shock  wave,  Chen  this 
[ shock  may  be  approximated  by  a shock  wave  normal  to  the  free  stream  flow.  To 

[ be  consistent  with  this  approximation  we  must  also  assume  that  the  motion  of 

i 
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any  shock  wave  that  arises  from  unsteady  changes  In  the  flow,  as  well  as  the 
motion  of  existing  shock  waves.  Is  also  calculated  by  this  normal  shock 
approximation.  For  this  simplified  model,  (3)  and  (4)  reduce  to 


1(^1  - 0 on  X = (dx/dt)  =»  ^ 

s s iK 


(y  + DM? 
00 


(10) 


which  gives  the  speed  of  the  normal  shock  In  the  flow  field.  For  steady  flows 
Is  a function  of  x alone;  this,  of  course,  still  permits  [(|i  1 to  vary 

X X 

• 

with  y.  For  unsteady  flows,  while  x Is  a function  of  t alone,  the 

s 

strength  of  the  shock  will  still  vary  with  y. 

For  tlme-llnearlzed  flows  the  steady  state  result  for  <|)®  with  normal 
shock  fitting  will  give  a steady  state  shock  position  x*,  0 i |y|  ^ y*  • 

We  now  determine  the  shock  wave's  motion  by  writing  the  perturbed  shock  posi- 
tion as  X ■ X®  + 6x(t)  and  using  the  time-linearized  version  of  (10);  we 
s s 

note  that  an  expression  of  the  form  x®  + 6x(t)/K  would  probably  be  more 
appropriate.  From  (10) , we  conclude  that  the  shock  motion  is  governed  by 


^ ■ ^2k  ^ 'l'jj(x,0,t)  with  1(^1  « I(|i®I  + fiE'I'l  ■ 0 


(11) 


on  the  shock.  Linearizing  the  expression  In  (11)  for  the  velocity  potential 
about  the  steady  shock  position  we  find 


(Kx  ,y,t)  - <|i(x®,y,t)  + ((>  (x®,y,t)<5dx 

9 S X S 

- 't'’’(x®,y)  + (j»®(x®,y)6dx  + <5'l'(x®,y,t)  + O(d^). 

9X9  5 


Because  we  have  treated  the  shock  as  a normal  one,  y appears  here  simply  as 
a parameter.  Now  IiKx^,y,t)l  and  Iij>®(x®,y)I  are  both  zero;  consequently 


we  have 


9 


ft.  . . 

* ® ;o  * ® 


which  must  be  integrated  in  time  in  conjunction  with  the  solution  to  (8) . 


3.  NUMERICAL  PROCEDURES 


In  a preliminary  study  of  the  unsteady  transonic  flows  a normal  shock- 


fitting procedure  was  implemented  in  an  implicit-iterative  scheme.  Satis- 


factory results  were  obtained,  but  the  procedure  was  time-consuming  because 


of  the  iterative  process  required  at  each  time  step.  The  recent  studies  of 


Ballhaus  and  Steger  (11)  and  Ballhaus  and  Goorjian  (22)  show  that  an  ADI 


scheme  is  more  efficient  than  the  implicit- iterative  scheme  in  treating  the 


low  frequency  transonic  flows.  The  shock-fitting  algorithm  was  modified 


and  Implemented  with  an  ADI  scheme.  In  this  section  the  ADI  procedure  and 


the  method  used  for  unsteady  shock  fitting  are  briefly  reviewed. 


3.1  Coordinate  Stretching 


To  minimize  the  far-fleld  boundary  effects  on  the  numerical  results  a 


relatively  large  computational  region  is  usually  required.  For  some  of  the 


cases  studied  in  this  paper  the  shock  excursions  are  large  and  the  unsteady 


disturbances  carried  several  chord  lengths  away  from  the  airfoil;  thus,  the 


use  of  a relatively  large  computational  domain  seems  desirable.  A simple 


and  straightforward  way  of  computing  the  solution  in  a large  computational 


domain  is  to  use  nonuniform  mesh  distributions  with  most  of  the  mesh  points 


concentrated  in  the  region  of  interest.  An  alternative  is  to  Introduce 


analytical  coordinate  stretchings.  In  the  present  study,  we  use  the  following 


coordinate  stretchings: 


I 
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C ■ ±{l  - exp(+aj^x)}  for  x < 0 and  rj  ■ ±{1  - expC+a^y))  for  y < 0, 


where  a^  and  a^  are  constanta  that  control  the  mesh  distributions.  The 
Infinite  physical  domain  Is  transformed  Into  the  finite  computational  domain 
bounded  by  |S|  ^ 1>  and  |n|  ± 1.  The  transformation  provides  a concen- 
trated mesh  distribution  near  the  airfoil  which  Is  suitable  for  the  present 
study.  While  this  scaling  Is  not  consistent  with  Che  known  algebraic  decay 
of  Che  perturbations,  calculations  made  with  an  algebraic  scaling,  viz., 

£ - x/(lxl  + a^)  etc.,  gave  essentially  Identical  results.  The  exponential 
variation  used  here  seems  more  desirable  near  the  airfoil. 

The  governing  equation  (1),  written  In  the  stretched  coordinate  system. 
Is 


-2KM 


a^l-  Inl)  , 


- 


(Y  + 1)M^ 

2&hi  - |n|) 


^ M^-l 


(Y  + DM' 


a^d  - )5|)^^ 


. 1 - Ini  , 

\ i^ri  - TCI)'  ^ 


- 0. 


(13) 


Because  (13)  Is  In  divergence-free  form,  a conservative  difference  approxima- 
tion can  be  constructed  If  the  shock  wave  Is  to  be  "captured"  rather  Chan 
"fitted." 


The  normal  shock  Jump  relation  follows  directly  from  (13);  this  relation 
and  the  boundary  condition  on  the  airfoil  surface  are  now 


(d£/dt)^  - 


a^(l  - |£|)(y  + 1) 


2K 


- 1 


(y  + DM. 


j + a^(l  - |£|)* 


(lA) 


and 
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n a-  34  a-  at 


-1  + e-rp(a. /2)  ^ C 1 - exp(-a. /2). 


Equations  analogous  to  (13)  and  (14)  for  the  time-linearized  results  are 
given  in  Reference  10. 


3. 2 Alternating-Direction  Implicit  (ADI)  Method 

The  low  frequency  equation  In  the  stretched  coordinate  system  is 
solved  by  the  alternating-direction  implicit  scheme  developed  by  Ballhaus 
and  Steger  (11).  To  simplify  this  discussion,  equation  (13)  is  rewritten  in 
the  form 


^t  ^ ^ 


where  the  function  V,  F and  G may  be  determined  by  comparing  equations 
(13)  and  (15).  The  solution  is  advanced  from  time  level  "n"  to  level 
"n+1"  by  the  following  two-step  procedure: 


^ - H'“)  + D^f"^  + - 0;  ^ ('T^^  ‘ I 


Hare  'V  refers  to  an  intermediate  value  of  V,  D^  is  the  type-dependent 
difference  operator  for  4~<lerlvatlves  and  the  central-difference  approx- 
imation for  n-derlvatlve.  The  baclcward  difference  approximation  for  can 

be  either  a first-order  or  a second-order  difference  approximation,  with  the 
latter  giving  Improved  results.  The  nonlinear  term  F is  evaluated,  using  a 
linearization  somewhat  different  from  the  two-time  level  averaging  procedure 


ii«i’rivi«0ii  «wi 


of  Ballhaus  and  Steger.  The  difference  approximations  described  above  pro- 


vide first-  or  second-order  accuracy  for  second-order  accuracy  for 

and  In  subsonic  regions,  and  first-order  accuracy  for  in 

supersonic  regions.  A local  analysis  shows  that  the  procedure  Is  uncondi- 
tionally stable. 

In  the  first  step  a quadradlagonal  system  is  generated  and  can  be  easily 


solved  by  direct  elimination.  For  lifting  calculations  two  grid  lines  are 
used  to  represent  the  lower  and  upper  surfaces  of  the  airfoil.  The  circula- 
tion, r.  Is  calculated  by  T ■ ~ ♦tttt  through  each  sweep.  Here 

"ITE"  denotes  the  upper  and  lower  values  at  the  first  grid  point  behind  the 
trailing  edge.  This  circulation  Is  Incorporated  Into  the  construction  of 


the  n-derivatlves  behind  the  airfoil  for  n ••  0, 


In  the  second  step  a tridiagonal  system  Is  generated  by  the  body.  Ahead 
of  the  leading  edge  and  behind  the  trailing  edge  the  double  grid  notation  for 
n ■ 0 destroys  the  trldlagonal  system.  However,  ahead  of  the  leading  edge, 
and  behind  the  trailing  edge,  + F;  thus,  the  difference 


equations  can  be  reordered  to  give  a trldlagonal  system.  On  the  airfoil  sur- 
face, the  matrix  equations  above  and  below  the  airfoil  are  decoupled;  they 
can  either  be  solved  separately  or  simultaneously  by  packing  the  matrix  equa- 
tions together. 


Again,  analogous  but  somewhat  simpler  equations  and  procedures  are  used 
for  the  time- linearized  calculations.  In  these  calculations  the  type  depen- 
dent operator,  D^,  changes  at  the  steady  state  sonic  line  and  shock  wave. 
The  coefficient,  f(C,n),  that  appears  in  (8)  in  the  form  {f 
depends  on  the  steady  state  results  end  must  be  stored.  On  the 

other  hand,  the  matrices  used  do  not  depend  on  the  solution  ip  and,  conse- 
quently, need  only  be  Inverted  once.  In  Its  present  form  our  algorithm  does 
not  take  advantage  of  this  feature. 


t 

[ 

I 

\ 

I 
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3.3  Shock  Fitting 

We  start  the  unsteady  nonlinear  flow  calculations  by  using  an  ADI 
scheme.  When  the  local  flow  becomes  sonic  and  compressive,  we  Introduce 
the  shock-fitting  algorithm  described  in  detail  in  Reference  21.  Sonic, 
compressive  points  are  treated  as  shock  points  where  differentiation  in  t 
and  ^ across  discontinuities  is  avoided.  Initially,  the  shock  has  zero 
strength  and  is  stationary.  The  flow  properties  ahead  of  and  behind  the 
shock  can  be  easily  extrapolated  from  neighboring  points.  The  shock  wave 
can  either  Increase  or  decrease  in  strength  during  the  unsteady  process. 

This  results  in  three  possibilities  for  shock  motion  that  have  to  be  con- 
sidered separately  in  the  fitting  procedure:  The  shock  moves  upstream  and 
crosses  grid  points;  the  shock  remains  stationary  or  moves  within  a grid 
spacing;  the  shock  moves  downstream  and  crosses  grid  points.  At  each  new 
time  level  the  shock  position  is  determined  by  applying  (10) . The  formula- 
tions of  the  difference  approximations  for  each  case  are  quite  similar. 

For  time-linearized  calculations  the  solution  is  advanced  in  time  using 
the  time-linearized  analogues  of  (16)  coupled  with  (12)  in  the  form 

Nl'*’  - -C(n)At((»“  + Hi”;  ■ -c(n)Ati(»^^  + Hi"*". 

Here 

c(n)  |c;|)^I^J(5;,n)I, 

and  C”  denotes  the  steady  state  position  of  the  shock  wave.  This  procedure 
corrects  the  4*  values  for  shock  motions  as  the  solution  progresses.  The 
shock  motion  is  easily  determined  simultaneously  by  using  (11)  and  (12)  in 
the  form 


I 
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0)1. 

9 A S 

Further  details  are  given  In  Reference  10. 


4.  RESULTS  AND  DISCUSSION 

Both  the  nonlinear  and  time- linearized  algorithms  have  been  used  to 
compute  the  flow  past  an  NACA  64A006  airfoil  subjected  to  Indlclal,  l.e., 
step,  changes  and  harmonic  motions  In  pitch  and  flap  oscillation.  The  lat- 
ter calculations  have  included  a range  of  Mach  numbers,  amplitudes  for  the 
nonlinear  algorithm,  and  reduced  frequencies  for  the  harmonic  changes.  The 
nonlinear  algorithm  has  also  been  used  to  compute  the  flow  past  a pulsating 
parabolic  arc  airfoil.  In  this  latter  flow,  at  • 0.85,  as  the  airfoil 
thickens  a shock  wave  forms  and  moves  downstream  until  shortly  after  mid- 
cycle. As  the  airfoil  thins,  the  shock  wave  moves  upstream  with  Increasing 
speed,  eventually  leaving  the  airfoil.  A comparison  of  the  results,  with 
and  without  shock  fitting  (21),  Indicates  that  shock  fitting  predicts  the 
formation  of  the  shock  wave  more  accurately.  It  also  properly  defines  the 
shock  wave  when  It  becomes  "subsonic-subsonic"  In  the  fixed  grid  system.  The 
shock  wave  decays  slowly  as  It  propagates  Into  the  free  stream  after  passing 
the  location  of  the  leading  edge  when  the  airfoil's  thickness  has  Just  become 
zero. 


4. 1 NACA  64A006  Airfoil.  Nonlinear  Calculations 

Steady  state  solutions  were  computed  as  discussed  In  Section  2.1  for  an 
NACA  64A006  airfoil  for  various  values  of  the  freestream  Mach  number  by  usln^ 
the  ADI  scheme  with  shock  fitting  outlined  In  Section  3.  The  free  stream 
Mach  number  was  varied  between  0.8  and  0.9.  The  mesh  system  had  101  by  82 
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grld  points  In  the  x-  and  y-dlrectlons  respectively.  About  250  to  450 
time  steps  were  required  for  the  solution  to  converge  | A(^  ^ 10 

These  steady  state  solutions  are  used  as  Initial  data  for  the  nonlinear  and 
tlme-llnearlzed  unsteady  flow  calculations. 

Results  were  computed  for  the  airfoil  with  quarter-chord  flap  for 
various  values  of  the  reduced  frequency,  the  free  stream  Mach  number,  and 
the  oscillation  amplitude.  In  order  to  simulate  the  shock  motions  observed 
by  Tljdeman  (5,  6).  These  motions  were  classified  by  him  as:  type  A - small 
shock  oscillation;  type  B - the  shock  becomes  very  weak  or  disappears  during 
part  of  a cycle;  type  C - the  shock  leaves  the  airfoil.  Results  for  type  A 
motions  are  not  given,  as  they  are  easy  to  treat  computationally.  For  all 
cases  studied  It  took  three  to  six  cycles  for  the  flow  field  to  become  peri- 
odic. Stability  seems  to  require  that  the  time  step  be  small  enough  that 
At (In  degrees) /K  < 10. 

Figure  1 Illustrates  the  pressure  coefficients  on  the  airfoil  surface  | 

j 

at  various  times  for  • 0.854,  K ■ 0.358  and  6 « 1“.  For  these  condl-  | 

I 

tlons  Ballhaus  and  Goorjlan  (12)  were  able  to  simulate  type  B motion  where  | 

the  shock  disappears  during  some  part  of  the  cycle.  Here  the  shock  does  not  ; 

disappear  during  the  cycle;  rather.  It  becomes  quite  weak  during  a small 

i 

portion  of  the  cycle.  This  difference  Is  probably  due  In  part  to  the  assump- 
tion of  a normal  shock,  which  results  In  a stronger  shock  than  would  normally  i 

i 

occur,  and  to  the  use  of  shock  fitting,  which  Is  able  to  resolve  very  weak  | 

shock  waves. 

Figure  2 depicts  the  pressure  coefficient  on  the  airfoil  surfaces  for 
* 0.822,  K ■ 0.496  and  d > 2",  simulating  type  C shock  motion.  Because 
we  have  used  less  spatial  resolution  and  have  not  scaled  the  equation  and 
boundary  conditions  with  various  powers  of  the  Mach  number,  a slightly  larger 
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deflectlon  angle  seems  to  be  needed  In  order  to  generate  the  type  C shock 
motion;  that  Is,  we  need  a 2*  deflection  angle  rather  than  the  1.5°  of 
Reference  12  to  obtain  analogous  behavior.  In  this  case  the  flow  field  Is 
subcrltlcal  during  most  of  the  cycle,  where  the  shock  wave  Is  barely  "cap- 
tured" In  the  non-shock-f Ittlng  procedure.  During  the  unsteady  process  the 
shock  moves  toward  the  leading  edge.  However,  the  strong  singular  behavior 
In  pressure  at  the  leading  edge  prevents  the  shock  from  propagating  off  the 
airfoil.  The  perturbation  velocity  becomes  large  and  Is  negative;  thus,  the 
flow  used  to  calculate  the  relative  velocity  ahead  of  the  shock  can  no  longer 
support  a shock  wave.  Normal  shock-fitting  calculations  determine  the  shock 
speed  from  the  pressure  jump  across  the  shock  at  the  airfoil  surface.  This 
eliminates  the  possibility  that  a portion  of  the  shock  may  propagate  off  the 
leading  edge  In  the  computations.  But  this  does  not  Imply  It  cannot  occur; 
rather  this  Is  a limitation  of  Che  normal  shock  fitting. 

Magnus  and  Yoshlhara  (15)  have  solved  the  Euler  equations  using  an 
explicit  procedure  for  Che  conditions  of  Figure  1.  Their  results  are  compared 
with  our  calculation  In  Figure  3 for  two  angular  times  chosen  to  represent  the 
least  and  Che  largest  discrepancies.  These  discrepancies  are  thought  Co  be 
mainly  due  to  Che  Inaccuracy  of  the  small  perturbation  solution  near  the 
leading  edge.  Small  errors  there  change  the  size  and  shape  of  Che  sonic  line 
and  influence  the  shock's  position.  For  the  conditions  considered,  the  shock 
Is  nearly  normal  and  the  normal  shock  approximation  should  be  a good  one. 
Rather  good  agreement  Is  obtained. 

Additional  nonlinear  calculations  have  been  carried  out  for  M > 0.880 
and  K 0.48.  Both  pitching  and  flap  motions  have  been  calculated  for  Indl- 
clal  and  harmonic  changes.  For  these  conditions  very  small  unsteady  changes 
lead  CO  very  small  shock  motions  and  Che  shock  wave  remains  becweed  grid 
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polnts.  Because  of  the  extrapolation  procedure  used  In  the  shock-fitting, 
the  ^ mesh  distribution  used  here  can  Introduce  errors,  albeit  small  ones, 
in  the  shock's  position  when  a grid  line  is  crossed.  We  wished  to  eliminate 
these  errors  in  order  to  use  the  nonlinear  calculations  to  judge  the  accuracy 
of  time-linearized  calculations.  These  results  indicate  that  for  pitching 
about  mid-chord,  nonlinear,  amplitude  dependent,  behavior  occurs  for 
6/t  ^ 0.1  for  K ■ 0.48.  Because  the  amplitude  of  the  shock  motions  in- 
creases with  decreasing  K,  nonlinear  effects  occur  at  smaller  values  of 
6/t  at  lower  reduced  frequencies. 

Indiclal  motions  require  about  eight  hundred  time  steps  of  varying  size 
to  resolve  the  response.  For  harmonic  motions,  initiated  from  rest,  three  to 
ten  cycles  are  required  for  the  solution  to  become  harmonic,  with  large  values 
of  K and  requiring  more  cycles.  The  pitching  mode  requires  more  cycles 

than  the  flap  mode.  The  amplitude  of  the  positive  and  negative  phases  of  the 
motion  could  be  varied  from  cycle  to  cycle  to  reduce  the  number  of  cycles 
required.  Each  cycle  requires  60  to  180  time  steps  to  compute,  with  more 
steps  required  for  smaller  values  of  K.  Each  time  step  takes  about  5 seconds 
of  CPU  time  on  a CDC  6400,  or  about  0.25  seconds  of  CPU  time  on  a CDC  7600. 


4.2  SAGA  64A006  Airfoil.  Time-Linearized  Calculations 

Time-linearized  results  have  been  computed  for  an  NACA  64A006  airfoil 
experiencing  harmonic  and  indiclal  pitching  and  flap  motions.  A noted  earlier, 
in  the  low  frequency  approximation  made  here  pitching  and  plunging  motions 
lead  to  the  same  result  except  that  the  time-linearized  perturbations  are 
proportional  to  the  maximum  pitch  angle  for  the  former,  and  K times  the 
maximum  amplitude  for  the  latter.  Harmonic  motions  initiated  from  a steady 
state  become  nearly  periodic  in  three  to  ten  cycles,  with  the  changes  induced 


by  flap  oscillations  becoming  periodic  more  rapidly  than  those  resulting 
from  pitching  oscillations.  More  cycles  were  required  for  larger  reduced 
frequencies  and,  to  a lesser  degree,  higher  Mach  numbers. 

In  order  to  confirm  the  validity  of  the  tlme-llnearlzed  calculations, 
both  the  tlme-llnearlzed  and  nonlinear  algorithms  were  used  to  compute  the 
response  to  a step  change  In  angle  of  attack  and  the  harmonic  response  to 
pitching  motions.  Figure  4 compares  the  nonlinear  and  tlme-llnearlzed  re- 
sults for  the  normalized  circulation  and  shock  position  for  harmonic 
pitching  motions  at  M^  ■>  0.88  and  K > 0.48.  Results  are  given  for  the 
fifth  cycle;  note  that  the  nonlinear  results  are  not  yet  periodic.  Figure 
5 compares  the  nonlinear  and  tlme-llnearlzed  pressure  deviation  from  steady 
state  at  six  angular  times  for  the  same  conditions.  Good  agreement  between 
Che  results  Is  obtained  for  6/t  less  than  0.1. 

Tlme-llnearlzed  pressure  distributions  at  six  angular  positions  for  an 
oscillating  quarter-chord  flap  with  K - 0.06  and  M^  » 0.875  are  shown  in 
Figure  6.  The  flap  deflection  Is  downward  during  Che  first  half  of  the 
cycle.  The  results  for  Che  second  half  of  the  period,  for  Che  symmetrical 
problem  shown  here,  are  Just  the  results  shown  with  Che  lower  and  upper 
surface  pressures  Interchanged.  Thus  Che  results  for  0°  are  not  given  as 
they  are  Just  above  those  for  180"  with  the  lower  and  upper  surface  pres- 
sures reversed.  Because  the  flap  hinge  occurs  very  close  to  the  steady 
state  shock  location,  the  pressure  singularity  due  to  Che  change  In  flow 
direction  at  the  hinge  Is  missed.  The  circulation  and  shock  excursion  obey 
Che  following  relations: 

r(t)/6  - 9.26  sin  (t  - 59"), 


X(t)  - 12  sin  (t  - 51"). 
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Note  the  substantial  phase  lag  In  the  circulation  and  the  shock's  position. 

Tlme-llnearized  pressure  distributions  at  six  angular  positions  for  an 
oscillating  airfoil  with  K * 0.12  and  « 0.875  are  depicted  in  Figure 
7.  If  these  results  are  multiplied  by  K,  then  they  represent  the  pressure 
perturbations  for  a plunging  airfoil.  As  in  the  previous  case  of  an  oscil- 
lating flap,  changes  In  forces  and  moments  of  0(5/K)  occur  due  to  shock 
wave  motion.  In  this  case 

r(t)/6  - 5.48  sin  (t  - 70"), 
x(t)  ■ 5.62  sin  (t  - 87"). 

Analogous  computations  have  been  carried  out  for  K >■  0.12,  0.24,  0.36, 
and  0.48.  Figure  8 depicts  the  shock  wave's  exucrslon  and  maximum  circula- 
tion as  a function  of  K The  nearly  linear  variation  of  the  shock 

excursion  substantiates  an  observation  made  In  a one-dimensional  model  where 
the  shock  wave  excursion  Is  directly  proportional  to  1/K  (see  Reference  10). 

In  these  calculations  the  circulation  gives  an  Immediate  evaluation  of 

the  lift  coefficient  as  a function  of  time;  the  moment  coefficient  must  be 

evaluated  by  Integrating  the  moment  of  the  pressure  coefficient.  This  is 

done  by  Integrating  the  moment  of  pressure  perturbations  with  the  shock  wave 

In  Its  steady-state  position  and  then  correcting  these  results  for  the  moment 

due  to  the  shock  wave  motion,  assuming  that  the  shock's  strength  Is  defined 

by  the  steady-state  pressure  field.  This  makes  an  error  In  the  shock  strength 

2 

of  0(6),  but  the  effect  on  the  moment  Is  0(6  /K) ; because  we  have  neglec- 
ted other  higher-order  terms  It  Is  consistent  to  neglect  this  change  In  the 
strength  of  the  shock  wave. 
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Flgure  9 depicts  the  absolute  value  and  phase  angle  of  the  normalized 
lift  and  moment  coefficients,  as  a function  of  the  Inverse  reduced  frequency 
K for  harmonic  flap  and  pitching  motions  at  ••  0.875. 

The  tlme-llnearlzed  algorithm  used  here  Is  a derivative  of  that  used 
for  the  nonlinear  calculations.  Consequently,  computational  times  are  not 
greatly  reduced  from  those  required  for  the  nonlinear  calculations.  The 
linearity  of  these  computations  may  make  It  possible  to  greatly  reduce  the 
computational  effort  required.  A local  stability  analysis  shows  that  the 
computations  should  be  unconditionally  stable,  but  numerical  experience  has 
shown  some  difficulties  for  At (In  degrees) /K  50.  Each  time  step  requires 
about  two  seconds  of  CPU  time  on  a CDC  6400,  or  about  0.1  seconds  on  a CDC 
7600.  The  number  of  time  steps  required  for  a given  computation  Is  somewhat 
less  than  those  required  for  the  nonlinear  computations  at  small  values  of 
K,  and  comparable  at  larger  values  of  K. 


5.  CONCLUSION 

Efficient  and  accurate  methods  for  computing  low  frequency,  unsteady 
behavior  In  transonic  flows  have  been  developed.  They  utilize  the  ADI  pro- 
cedure developed  at  NASA  Ames  for  the  small  perturbation  equation,  but  treat 
shock  waves  as  discontinuities.  The  tlme-llnearlzed  calculations  allow 
shock  wave  motions,  which  are  shown  to  be  0(6/K)  and  often  dominate  changes 
In  the  force  and  moment  coefficients.  Comparison  of  the  tlme-llnearlzed 
results  with  fully  nonlinear  calculations  delineates  their  range  of  appli- 
cability. The  unsteady  behavior  due  to  harmonic  pitching  and  flap  oscilla- 
tions of  an  NACA  64A006  airfoil  Is  discussed. 
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FICUU  S.  NOSMALIZED  MORLINEAk  (*•**)  AND  TIMB-LINKAXIZED  ( ) PUSSUKE 

PEXTUmilONS  ON  THE  (JPPEH  SUHEACE  OF  AN  NACA  6AA006  AT  SIX  TIMES 
FITCHINC  MOTION  WITH  M.  - 0.880,  t - 0.A8.  FOR  THE  NONLINEAR 
CALCULATIONS  < « 0.1*. 


